Introduction

…

Packages

library(lidR)# for working with  LAZ files
library(sf) # for working spatial class
library(raster)# for working with raster
library(rayshader) # for 3D viz
library(rgl) # for interactive plots

Motte

motte_point = st_read("data/canmore_point.shp", quiet = TRUE)

motte_buffer = st_buffer(motte_point,dist = 50)

Point Cloud

Next, I am reading the LAZ files and clipping the point cloud to the extent of immidiate area around the motte.

#read laz files
las = readLAS("data/NX6055_4PPM_LAS_PHASE3.laz")

motte = clip_roi(las, motte_buffer)
#plot lidar point cloud
plot(motte, bg = "white")

rglwidget()

grab and rotate the model !

Digital Elevation Model

In this step I am running standard algorithms from LiDR package to compute Digital Surface Model (DSM) and Digitial Terrain Model (DTM).

rgl.clear()
# create dsm and dtm
dsm = grid_canopy(motte, 0.5, pitfree())

# assign coordinate system
crs(dsm) = CRS("+init=epsg:27700")

# create dtm
dtm = grid_terrain(motte, 0.5, algorithm = knnidw(k = 6L, p = 2))

# addign coordinate system
crs(dtm) = CRS("+init=epsg:27700")

par(mfrow = c(1,2))
plot(dtm, main = "DTM", col = height.colors(50))
plot(dsm, main = "DSM",col = height.colors(50))

Hillshade

# dsm
slope_dsm = terrain(dsm, opt = 'slope')
aspect_dsm = terrain(dsm, opt = 'aspect')
hill_dsm = hillShade(slope_dsm, aspect_dsm, angle = 40, direction = 270)

# dtm
slope_dtm = terrain(dtm, opt = 'slope')
aspect_dtm = terrain(dtm, opt = 'aspect')
hill_dtm = hillShade(slope_dtm, aspect_dtm, angle = 5, direction = 270)

#plot
par(mfrow = c(1,2))
plot(hill_dtm, main = "DTM Hilshade", col = grey.colors(100, start = 0, end = 1), 
     legend = FALSE)
plot(hill_dsm, main = "DSM Hillshade", col = grey.colors(100, start = 0, end = 1))

3D Model

#And convert it to a matrix:
elmat = raster_to_matrix(dtm)

elmat %>%
  sphere_shade(texture = "imhof1") %>%
  add_shadow(ray_shade(elmat, zscale = 0.5, sunaltitude = 30,lambert = TRUE),
             max_darken = 1) %>%
  add_shadow(lamb_shade(elmat,zscale = 0.5,sunaltitude = 30), max_darken = 0.2) %>%
  add_shadow(ambient_shade(elmat, zscale = 0.5), max_darken = 0.2) %>%
  plot_3d(elmat, zscale = 0.5,windowsize = c(1000,1000),
          phi = 40, theta = 180, zoom = 0.8, fov = 1)
rglwidget()

grab and rotate the model !

Advanced Example

Tantallon Castle And here is the code

library(rayshader)
library(raster)
library(magick)
library(gifski)

dem = raster("data/tantallon.tif")

n_frames <- 41

elmat = raster_to_matrix(dem)

zscale <- 0.5

waterdepthvalues <- c(0:20, seq(19,0,-1)) / 2

phi_values = seq(20,60)

waterdepthvalues * zscale

length(waterdepthvalues)

img_frames <- paste0("drain", seq_len(n_frames), ".png")

for (i in seq_len(n_frames)) {
elmat %>%
  sphere_shade(texture = "imhof1") %>%
  add_shadow(ray_shade(elmat, zscale = 0.5, sunaltitude = 30,lambert = TRUE),
             max_darken = 1) %>%
  add_shadow(lamb_shade(elmat,zscale = 0.5,sunaltitude = 30), max_darken = 0.2) %>%
  add_shadow(ambient_shade(elmat, zscale = 0.5), max_darken = 0.2) %>%
  plot_3d(elmat, zscale = 0.5,windowsize = c(1000,1000),
          water = TRUE, watercolor = "imhof3", wateralpha = 0.5, 
          waterlinecolor = "#ffffff", waterlinealpha = 0.5,
          waterdepth = waterdepthvalues[i],
          phi = 50, theta = 45, zoom = 0.8, fov = 1)
  render_snapshot(filename = img_frames[i],
                  title_text = "Tantallon Castle", 
                  title_font = "Helvetica",
                  vignette = TRUE,
                  title_size = 50,
                  title_color = "black")
  rgl::clear3d()
}

magick::image_write_gif(magick::image_read(img_frames), 
                        path = "tantallon.gif", 
                        delay = 6/n_frames)

References

Roussel, Jean-Romain, and David Auty. 2021. Airborne LiDAR Data Manipulation and Visualization for Forestry Applications. https://cran.r-project.org/package=lidR.
Roussel, Jean-Romain, David Auty, Nicholas C. Coops, Piotr Tompalski, Tristan R. H. Goodbody, Andrew Sánchez Meador, Jean-François Bourdon, Florian de Boissieu, and Alexis Achim. 2020. “lidR: An r Package for Analysis of Airborne Laser Scanning (ALS) Data.” Remote Sensing of Environment 251: 112061. https://doi.org/10.1016/j.rse.2020.112061.